FEM-BEM Coupling

The ngbem boundary element addon project initiated by Lucy Weggeler (see https://weggler.github.io/ngbem/intro.html) is now partly integrated into core NGSolve. Find a short and sweet introduction to the boundary element method there.

In this demo we simulate a plate capacitor on an unbounded domain.

from ngsolve import *
from netgen.occ import *
from ngsolve.solvers import GMRes, MinRes
from ngsolve.webgui import Draw
from ngsolve.bem import *

largebox = Box ((-2,-2,-2), (2,2,2) )
eltop = Box ( (-1,-1,0.5), (1,1,1) )
elbot = Box ( (-1,-1,-1), (1,1,-0.5))

largebox.faces.name = "outer" # coupling boundary
eltop.faces.name = "el1" # Dirichlet boundary
elbot.faces.name = "el2" # Dirichlet boundary
eltop.edges.hpref = 1
elbot.edges.hpref = 1

shell = largebox-eltop-elbot # FEM domain 
shell.solids.name = "air"

mesh = shell.GenerateMesh(maxh=0.8)
mesh.RefineHP(2)
ea = { "euler_angles" : (-67, 0, 110) } 
Draw (mesh, clipping={"x":1, "y":0, "z":0, "dist" : 1.1}, **ea);

On the exterior domain \(\Omega^c\), the solution can be expressed by the representation formula:

\[ x \in \Omega^c: \quad u(x) = - \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{1}{\| x-y\|} } \, \gamma_1 u (y)\, \mathrm{d}\sigma_y + \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{\langle n(y) , x-y\rangle }{\| x-y\|^3} } \, \gamma_0 u (y)\, \mathrm{d}\sigma_y\,, \]

where \(\gamma_0 u = u\) and \(\gamma_1 u = \frac{\partial u}{\partial n}\) are Dirichlet and Neumann traces. These traces are related by the Calderon projector

\[ \begin{align}\begin{aligned}\begin{split} \left( \begin{array}{c} \gamma_0 u \\ \gamma_1 u \end{array}\right) = \left( \begin{array}{cc} -V & \frac12 + K \\ \frac12 - K^\intercal & -D \end{array} \right) \left( \begin{array}{c} \gamma_1 u \\ \gamma_0 u \end{array}\right) $$.\end{split}\\The $V$, $K$ are the single layer and double layer potential operators, and $D$ is the hypersingular operator.\\On the FEM domain we have the variational formulation\end{aligned}\end{align} \]

\int_{\Omega_\text{FEM}} \nabla u \nabla v , dx - \int_\Gamma \gamma_1 u v , ds = 0 \qquad \forall , v \in H^1(\Omega_\text{FEM}) $$

We use Calderon’s represenataion formula for the Neumann trace:

\[ \int_{\Omega_\text{FEM}} \nabla u \nabla v \, dx - \int_\Gamma \left( \left( \tfrac{1}{2} - K^\intercal\right) \,\gamma_1 u - D \, \gamma_0 u\right) v = 0 \qquad \forall \, v \in H^1(\Omega_\text{FEM}) \]

To get a closed system, we use also the first equation of the Calderon equations. To see the structure of the discretized system, the dofs are split into degrees of freedom inside \(\Omega\), and those on the boundary \(\Gamma\). The FEM matrix \(A\) is split accordingly. We see, the coupled system is symmetric, but indefinite:

\[\begin{split} \left( \begin{array}{ccc } A_{\Omega\Omega} & A_{\Omega\Gamma} & 0 \\ A_{\Gamma\Omega} & A_{\Gamma\Gamma } + D & -\frac12 M^\intercal + K^\intercal \\ 0 & -\frac12 M + K & -V \end{array}\right) \left( \begin{array}{c} u \\ \gamma_0 u \\ \gamma_1 u \end{array}\right) = \left( \begin{array}{c} F_{\Omega} \\ F_{\Gamma}\\ 0 \end{array}\right) \,. \end{split}\]

Generate the finite element space for \(H^1(\Omega)\) and set the given Dirichlet boundary conditions on the surfaces of the plates:

order = 4
fesH1 = H1(mesh, order=order, dirichlet="el1|el2") 
print ("H1-ndof = ", fesH1.ndof)

H1-ndof =  91269

The finite element space \(\verb-fesH1-\) provides \(H^{\frac12}(\Gamma)\) conforming element to discretize the Dirichlet trace on the coupling boundary \(\Gamma\). However we still need \(H^{-\frac12}(\Gamma)\) conforming elements to discretize the Neumann trace of \(u\) on the coupling boundary. Here it is:

fesL2 = SurfaceL2(mesh, order=order-1, dual_mapping=True, definedon=mesh.Boundaries("outer")) 
print ("L2-ndof = ", fesL2.ndof)

L2-ndof =  4035

fes = fesH1 * fesL2
u,dudn = fes.TrialFunction()
v,dvdn = fes.TestFunction()

a = BilinearForm(grad(u)*grad(v)*dx, check_unused=False).Assemble() 

gfudir = GridFunction(fes)
gfudir.components[0].Set ( mesh.BoundaryCF( { "el1" : 1, "el2" : -1 }), BND)

f = LinearForm(fes).Assemble()
res = (f.vec - a.mat * gfudir.vec).Evaluate()

Generate the the single layer potential \(V\), double layer potential \(K\) and hypersingular operator \(D\):

n = specialcf.normal(3)
with TaskManager():
    V = LaplaceSL(dudn*ds("outer"))*dvdn*ds("outer")
    K = LaplaceDL(u*ds("outer"))*dvdn*ds("outer")
    D = LaplaceSL(Cross(grad(u).Trace(),n)*ds("outer"))*Cross(grad(v).Trace(),n)*ds("outer")
    M = BilinearForm(u*dvdn*ds("outer"), check_unused=False).Assemble()

Setup the coupled system matrix and the right hand side:

sym = a.mat+D.mat - (0.5*M.mat+K.mat).T - (0.5*M.mat+K.mat) - V.mat
rhs = res

bfpre = BilinearForm(grad(u)*grad(v)*dx+1e-10*u*v*dx  + dudn*dvdn*ds("outer") ).Assemble()
pre = bfpre.mat.Inverse(freedofs=fes.FreeDofs(), inverse="sparsecholesky")

Compute the solution of the coupled system:

with TaskManager():
    sol_sym = MinRes(mat=sym, rhs=rhs, pre=pre, tol=1e-8, maxsteps=200, printrates=True)

LinearSolver iteration 1, residual = 18.46865730865644     
LinearSolver iteration 2, residual = 2.3400628305255675     
LinearSolver iteration 3, residual = 2.082508582473063     

LinearSolver iteration 4, residual = 0.4214702285936898     
LinearSolver iteration 5, residual = 0.26625863646450654     

LinearSolver iteration 6, residual = 0.11205519209281295     
LinearSolver iteration 7, residual = 0.05786199306545903     

LinearSolver iteration 8, residual = 0.03529379998646658     
LinearSolver iteration 9, residual = 0.015622860652925912     

LinearSolver iteration 10, residual = 0.014077913354084292     
LinearSolver iteration 11, residual = 0.014053466804286464     
LinearSolver iteration 12, residual = 0.003550247631373225     

LinearSolver iteration 13, residual = 0.0030969657773826376     
LinearSolver iteration 14, residual = 0.002083719255311875     

LinearSolver iteration 15, residual = 0.0016137597476426266     
LinearSolver iteration 16, residual = 0.0013415624945746854     

LinearSolver iteration 17, residual = 0.0013275546909253192     
LinearSolver iteration 18, residual = 0.0007480877276348744     

LinearSolver iteration 19, residual = 0.0007138753088691281     
LinearSolver iteration 20, residual = 0.000570198420049705     

LinearSolver iteration 21, residual = 0.0005637143009037972     
LinearSolver iteration 22, residual = 0.0003837960885996448     

LinearSolver iteration 23, residual = 0.00038347012601368875     
LinearSolver iteration 24, residual = 0.0002819114975842227     

LinearSolver iteration 25, residual = 0.00027352917319117684     
LinearSolver iteration 26, residual = 0.00021513653437461877     

LinearSolver iteration 27, residual = 0.00020073809279816537     
LinearSolver iteration 28, residual = 0.000137966409497975     

LinearSolver iteration 29, residual = 0.00013768299243130408     
LinearSolver iteration 30, residual = 0.00012073734590311811     

LinearSolver iteration 31, residual = 0.0001076635832151948     
LinearSolver iteration 32, residual = 8.128459071219661e-05     

LinearSolver iteration 33, residual = 7.76683695611141e-05     
LinearSolver iteration 34, residual = 7.732407896670565e-05     

LinearSolver iteration 35, residual = 6.19761867588294e-05     
LinearSolver iteration 36, residual = 6.111656645258829e-05     

LinearSolver iteration 37, residual = 4.4089014379919995e-05     
LinearSolver iteration 38, residual = 4.003433127760303e-05     

LinearSolver iteration 39, residual = 3.706923103118659e-05     
LinearSolver iteration 40, residual = 3.062696319972191e-05     

LinearSolver iteration 41, residual = 2.7356170355283733e-05     
LinearSolver iteration 42, residual = 2.5400647581731705e-05     

LinearSolver iteration 43, residual = 2.2683596195650594e-05     
LinearSolver iteration 44, residual = 2.2566300870180387e-05     

LinearSolver iteration 45, residual = 2.033098562876504e-05     
LinearSolver iteration 46, residual = 1.931368872976361e-05     

LinearSolver iteration 47, residual = 1.7512760318945215e-05     
LinearSolver iteration 48, residual = 1.5749944383765617e-05     

LinearSolver iteration 49, residual = 1.2285656112928512e-05     
LinearSolver iteration 50, residual = 1.2239170472469926e-05     

LinearSolver iteration 51, residual = 1.0104200135183579e-05     
LinearSolver iteration 52, residual = 1.0043028891719694e-05     
LinearSolver iteration 53, residual = 9.837165878960726e-06     

LinearSolver iteration 54, residual = 8.01526970143426e-06     
LinearSolver iteration 55, residual = 7.837208847013061e-06     

LinearSolver iteration 56, residual = 6.749304895335784e-06     
LinearSolver iteration 57, residual = 6.5449044177439576e-06     

LinearSolver iteration 58, residual = 5.592815110851386e-06     
LinearSolver iteration 59, residual = 4.771397416191434e-06     

LinearSolver iteration 60, residual = 4.676201849053099e-06     
LinearSolver iteration 61, residual = 4.345523299371716e-06     

LinearSolver iteration 62, residual = 4.078718648048312e-06     
LinearSolver iteration 63, residual = 4.06216432838107e-06     

LinearSolver iteration 64, residual = 3.531416858269196e-06     
LinearSolver iteration 65, residual = 3.4010307839988133e-06     
LinearSolver iteration 66, residual = 2.7597114672827767e-06     

LinearSolver iteration 67, residual = 2.7596970962543564e-06     
LinearSolver iteration 68, residual = 2.386099685019108e-06     

LinearSolver iteration 69, residual = 2.3807107093287264e-06     
LinearSolver iteration 70, residual = 1.9815171659097434e-06     

LinearSolver iteration 71, residual = 1.9756841553701163e-06     
LinearSolver iteration 72, residual = 1.791964613515088e-06     

LinearSolver iteration 73, residual = 1.634351964035279e-06     
LinearSolver iteration 74, residual = 1.575189128239365e-06     

LinearSolver iteration 75, residual = 1.4355477508144318e-06     
LinearSolver iteration 76, residual = 1.4353299817311271e-06     

LinearSolver iteration 77, residual = 1.1026082333568388e-06     
LinearSolver iteration 78, residual = 1.0249144538761416e-06     

LinearSolver iteration 79, residual = 9.36516369425444e-07     
LinearSolver iteration 80, residual = 8.351570428810832e-07     

LinearSolver iteration 81, residual = 8.06946431604624e-07     
LinearSolver iteration 82, residual = 7.777619452585362e-07     

LinearSolver iteration 83, residual = 6.240082101816985e-07     
LinearSolver iteration 84, residual = 6.192931727982841e-07     

LinearSolver iteration 85, residual = 5.085600611168224e-07     
LinearSolver iteration 86, residual = 4.91738666983721e-07     

LinearSolver iteration 87, residual = 4.330863220186696e-07     
LinearSolver iteration 88, residual = 3.995443846956006e-07     

LinearSolver iteration 89, residual = 3.481635850900664e-07     
LinearSolver iteration 90, residual = 3.145969948339135e-07     

LinearSolver iteration 91, residual = 2.991235200445901e-07     
LinearSolver iteration 92, residual = 2.96942381125406e-07     

LinearSolver iteration 93, residual = 2.8463992146428413e-07     
LinearSolver iteration 94, residual = 2.608483427241906e-07     

LinearSolver iteration 95, residual = 2.4806545370355175e-07     
LinearSolver iteration 96, residual = 2.2154325166374918e-07     

LinearSolver iteration 97, residual = 2.0374354114640737e-07     
LinearSolver iteration 98, residual = 1.775163496575714e-07     

gfu = GridFunction(fes)
gfu.vec[:] = sol_sym + gfudir.vec 
Draw(gfu.components[0], clipping={"x" : 1, "y":0, "z":0, "dist":0.0, "function" : True }, **ea, order=2); 

The Neumann data:

Draw (gfu.components[1], **ea, order=3);